home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / mdot.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  14.4 KB  |  392 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;; (c) Copyright 1982 Massachusetts Institute of Technology 
  9.  
  10. (in-package "MAXIMA")
  11. ;; Non-commutative product and exponentiation simplifier
  12. ;; Written:    July 1978 by CWH
  13.  
  14. ;; Flags to control simplification:
  15.  
  16. (macsyma-module mdot)
  17.  
  18. (DEFMVAR $DOTCONSTRULES T
  19.      "Causes a non-commutative product of a constant and
  20. another term to be simplified to a commutative product.  Turning on this
  21. flag effectively turns on DOT0SIMP, DOT0NSCSIMP, and DOT1SIMP as well.")
  22.  
  23. (DEFMVAR $DOT0SIMP T
  24.      "Causes a non-commutative product of zero and a scalar term to
  25. be simplified to a commutative product.")
  26.  
  27. (DEFMVAR $DOT0NSCSIMP T
  28.      "Causes a non-commutative product of zero and a nonscalar term
  29. to be simplified to a commutative product.")
  30.  
  31. (DEFMVAR $DOT1SIMP T
  32.     "Causes a non-commutative product of one and another term to be
  33. simplified to  a commutative product.")
  34.  
  35. (DEFMVAR $DOTSCRULES NIL
  36.      "Causes a non-commutative product of a scalar and another term to
  37. be simplified to a commutative product.  Scalars and constants are carried
  38. to the front of the expression.")
  39.  
  40. (DEFMVAR $DOTDISTRIB NIL
  41.      "Causes every non-commutative product to be expanded each time it
  42. is simplified, i.e.  A . (B + C) will simplify to A . B + A . C.")
  43.  
  44. (DEFMVAR $DOTEXPTSIMP T "Causes A . A to be simplified to A ^^ 2.")
  45.  
  46. (DEFMVAR $DOTASSOC T
  47.      "Causes a non-commutative product to be considered associative, so
  48. that A . (B . C) is simplified to A . B . C.  If this flag is off, dot is
  49. taken to be right associative, i.e.  A . B . C is simplified to A . (B . C).")
  50.  
  51. (DEFMVAR $DOALLMXOPS T
  52.      "Causes all operations relating to matrices (and lists) to be
  53. carried out.  For example, the product of two matrices will actually be
  54. computed rather than simply being returned.  Turning on this switch
  55. effectively turns on the following three.")
  56.  
  57. (DEFMVAR $DOMXMXOPS T "Causes matrix-matrix operations to be carried out.")
  58.  
  59. (DEFMVAR $DOSCMXOPS NIL "Causes scalar-matrix operations to be carried out.")
  60.  
  61. (DEFMVAR $DOMXNCTIMES NIL
  62.      "Causes non-commutative products of matrices to be carried out.")
  63.  
  64. (DEFMVAR $SCALARMATRIXP T
  65.      "Causes a square matrix of dimension one to be converted to a
  66. scalar, i.e. its only element.")
  67.  
  68. (DEFMVAR $DOTIDENT 1 "The value to be returned by X^^0.")
  69.  
  70. (DEFMVAR $ASSUMESCALAR T
  71.      "This governs whether unknown expressions 'exp' are assumed to behave
  72. like scalars for combinations of the form 'exp op matrix' where op is one of
  73. {+, *, ^, .}.  It has three settings:
  74.  
  75. FALSE -- such expressions behave like non-scalars.
  76. TRUE  -- such expressions behave like scalars only for the commutative
  77.      operators but not for non-commutative multiplication.
  78. ALL   -- such expressions will behave like scalars for all operators
  79.      listed above.
  80.  
  81. Note:  This switch is primarily for the benefit of old code.  If possible,
  82. you should declare your variables to be SCALAR or NONSCALAR so that there
  83. is no need to rely on the setting of this switch.")
  84.  
  85. ;; Specials defined elsewhere.
  86.  
  87. (DECLARE-TOP(SPECIAL $EXPOP $EXPON ; Controls behavior of EXPAND
  88.           SIGN        ; Something to do with BBSORT1
  89.           ERRORSW)
  90.      (FIXNUM $EXPOP $EXPON)
  91.      (*EXPR FIRSTN $IDENT POWERX MXORLISTP1 ONEP1
  92.         SCALAR-OR-CONSTANT-P EQTEST BBSORT1 OUTERMAP1 TIMEX))
  93.  
  94. (defun simpnct (exp vestigial simp-flag) 
  95.  vestigial ;ignored
  96.  (let ((check exp)
  97.        (first-factor (simpcheck (cadr exp) simp-flag))
  98.        (remainder (if (cdddr exp)
  99.               (ncmuln (cddr exp) simp-flag)
  100.               (simpcheck (caddr exp) simp-flag))))
  101.       (cond ((null (cdr exp)) $dotident)
  102.         ((null (cddr exp)) first-factor)
  103.  
  104. ;  This does (. sc m) --> (f* sc m)  and  (. (f* sc m1) m2) --> (f* sc (. m1 m2))
  105. ;  and (. m1 (f* sc m2)) --> (f* sc (. m1 m2)) where sc can be a scalar
  106. ;  or constant, and m1 and m2 are non-constant, non-scalar expressions.
  107.  
  108.         ((commutative-productp first-factor remainder)
  109.          (mul2 first-factor remainder))
  110.         ((product-with-inner-scalarp first-factor)
  111.          (let ((p-p (partition-product first-factor)))
  112.           (outer-constant (car p-p) (cdr p-p) remainder)))
  113.         ((product-with-inner-scalarp remainder)
  114.          (let ((p-p (partition-product remainder)))
  115.           (outer-constant (car p-p) first-factor (cdr p-p))))
  116.  
  117. ;  This code does distribution when flags are set and when called by
  118. ;  $EXPAND.  The way we recognize if we are called by $EXPAND is to look at
  119. ;  the value of $EXPOP, but this is a kludge since $EXPOP has nothing to do
  120. ;  with expanding (. A (f+ B C)) --> (f+ (. A B) (. A C)).  I think that
  121. ;  $EXPAND wants to have two flags:  one which says to convert
  122. ;  exponentiations to repeated products, and another which says to
  123. ;  distribute products over sums.
  124.  
  125.         ((and (mplusp first-factor) (or $dotdistrib (not (zerop $expop))))
  126.          (addn (mapcar #'(lambda (x) (ncmul x remainder))
  127.                (cdr first-factor))
  128.            t))
  129.         ((and (mplusp remainder) (or $dotdistrib (not (zerop $expop))))
  130.          (addn (mapcar #'(lambda (x) (ncmul first-factor x))
  131.                (cdr remainder))
  132.            t))
  133.  
  134. ;  This code carries out matrix operations when flags are set.
  135.  
  136.         ((matrix-matrix-productp first-factor remainder)
  137.          (timex first-factor remainder))
  138.         ((or (scalar-matrix-productp first-factor remainder)
  139.          (scalar-matrix-productp remainder first-factor))
  140.          (simplifya (outermap1 'mnctimes first-factor remainder) t))
  141.  
  142. ;  (. (^^ x n) (^^ x m)) --> (^^ x (f+ n m))
  143.  
  144.         ((and (simpnct-alike first-factor remainder) $dotexptsimp)
  145.          (simpnct-merge-factors first-factor remainder))
  146.  
  147. ;  (. (. x y) z) --> (. x y z)
  148.  
  149.         ((and (mnctimesp first-factor) $dotassoc)
  150.          (ncmuln (append (cdr first-factor)
  151.                  (if (mnctimesp remainder)
  152.                  (cdr remainder)
  153.                  (ncons remainder)))
  154.              t))
  155.  
  156. ;  (. (^^ (. x y) m) (^^ (. x y) n) z) --> (. (^^ (. x y) m+n) z)
  157. ;  (. (^^ (. x y) m) x y z) --> (. (^^ (. x y) m+1) z)
  158. ;  (. x y (^^ (. x y) m) z) --> (. (^^ (. x y) m+1) z)
  159. ;  (. x y x y z) --> (. (^^ (. x y) 2) z)
  160.  
  161.         ((and (mnctimesp remainder) $dotassoc $dotexptsimp)
  162.          (setq exp (simpnct-merge-product first-factor (cdr remainder)))
  163.          (if (and (mnctimesp exp) $dotassoc)
  164.          (simpnct-antisym-check (cdr exp) check)
  165.          (eqtest exp check)))
  166.  
  167. ;  (. x (. y z)) --> (. x y z)
  168.  
  169.         ((and (mnctimesp remainder) $dotassoc)
  170.          (simpnct-antisym-check (cons first-factor (cdr remainder)) check))
  171.  
  172.         (t (eqtest (list '(mnctimes) first-factor remainder) check)))))
  173.  
  174. ;  Predicate functions for simplifying a non-commutative product to a
  175. ;  commutative one.  SIMPNCT-CONSTANTP actually determines if a term is a
  176. ;  constant and is not a nonscalar, i.e. not declared nonscalar and not a
  177. ;  constant list or matrix.  The function CONSTANTP determines if its argument
  178. ;  is a number or a variable declared constant.
  179.  
  180. (defun commutative-productp (first-factor remainder)
  181.        (or (simpnct-sc-or-const-p first-factor)
  182.        (simpnct-sc-or-const-p remainder)
  183.        (simpnct-onep first-factor)
  184.        (simpnct-onep remainder)
  185.        (zero-productp first-factor remainder)
  186.        (zero-productp remainder first-factor)))
  187.  
  188. (defun simpnct-sc-or-const-p (term)
  189.   (or (simpnct-constantp term) (simpnct-assumescalarp term)))
  190.  
  191. (defun simpnct-constantp (term)
  192.        (and $dotconstrules
  193.         (or (mnump term)
  194.         (and ($constantp term) (not ($nonscalarp term))))))
  195.  
  196. (defun simpnct-assumescalarp (term)
  197.        (and $dotscrules (scalar-or-constant-p term (eq $assumescalar '$all))))
  198.  
  199. (defun simpnct-onep (term) (and $dot1simp (onep1 term)))
  200.     
  201. (defun zero-productp (one-term other-term)
  202.        (and (zerop1 one-term)
  203.         $dot0simp
  204.         (or $dot0nscsimp (not ($nonscalarp other-term)))))
  205.  
  206. ;  This function takes a form and determines if it is a product 
  207. ;  containing a constant or a declared scalar.  Note that in the
  208. ;  next three functions, the word "scalar" is used to refer to a constant
  209. ;  or a declared scalar.  This is a bad way of doing things since we have
  210. ;  to cdr down an expression twice: once to determine if a scalar is there
  211. ;  and once again to pull it out.
  212.  
  213. (defun product-with-inner-scalarp (product)
  214.        (and (mtimesp product)
  215.         (or $dotconstrules $dotscrules)
  216.         (do ((factor-list (cdr product) (cdr factor-list)))
  217.         ((null factor-list) nil)
  218.         (if (simpnct-sc-or-const-p (car factor-list))
  219.             (return t)))))
  220.  
  221. ;  This function takes a commutative product and separates it into a scalar
  222. ;  part and a non-scalar part.
  223.  
  224. (defun partition-product (product)
  225.        (do ((factor-list (cdr product) (cdr factor-list))
  226.         (scalar-list nil)
  227.         (nonscalar-list nil))
  228.        ((null factor-list) (cons (nreverse scalar-list)
  229.                      (muln (nreverse nonscalar-list) t)))
  230.        (if (simpnct-sc-or-const-p (car factor-list))
  231.            (push (car factor-list) scalar-list)
  232.            (push (car factor-list) nonscalar-list))))
  233.  
  234. ;  This function takes a list of constants and scalars, and two nonscalar
  235. ;  expressions and forms a non-commutative product of the nonscalar
  236. ;  expressions, and a commutative product of the constants and scalars and
  237. ;  the non-commutative product.
  238.  
  239. (defun outer-constant (constant nonscalar1 nonscalar2)
  240.        (muln (nconc constant (ncons (ncmul nonscalar1 nonscalar2))) t))
  241.  
  242. (defun simpnct-base (term) (if (mncexptp term) (cadr term) term))
  243.  
  244. (defun simpnct-power (term) (if (mncexptp term) (caddr term) 1))
  245.  
  246. (defun simpnct-alike (term1 term2)
  247.        (alike1 (simpnct-base term1) (simpnct-base term2)))
  248.  
  249. (defun simpnct-merge-factors (term1 term2)
  250.        (ncpower (simpnct-base term1)
  251.         (add2 (simpnct-power term1) (simpnct-power term2))))
  252.  
  253. (defun matrix-matrix-productp (term1 term2)
  254.        (and (or $doallmxops $domxmxops $domxnctimes)
  255.         (mxorlistp1 term1)
  256.         (mxorlistp1 term2)))
  257.  
  258. (defun scalar-matrix-productp (term1 term2)
  259.        (and (or $doallmxops $doscmxops)
  260.         (mxorlistp1 term1)
  261.         (scalar-or-constant-p term2 (eq $assumescalar '$all))))
  262.  
  263. (declare-top(muzzled t))
  264.  
  265. (defun simpncexpt (exp vestigial simp-flag)
  266.   vestigial ;ignored
  267.   (let ((factor (simpcheck (cadr exp) simp-flag))
  268.     (power (simpcheck (caddr exp) simp-flag))
  269.     (check exp))
  270.        (twoargcheck exp)
  271.        (cond ((zerop1 power)
  272.           (if (zerop1 factor)
  273.           (if (not errorsw)
  274.               (merror "~M has been generated"
  275.                   (list '(mncexpt) factor power))
  276.               (throw 'errorsw t)))
  277.           (if (mxorlistp1 factor) (identitymx factor) $dotident))
  278.          ((onep1 power) factor)
  279.          ((simpnct-sc-or-const-p factor) (power factor power))
  280.          ((and (zerop1 factor) $dot0simp) factor)
  281.          ((and (onep1 factor) $dot1simp) factor)
  282.          ((and (or $doallmxops $domxmxops)
  283.            (mxorlistp1 factor)
  284.            (fixnump power) 'fixnum)
  285.           (let (($scalarmatrixp (or ($listp factor) $scalarmatrixp)))
  286.         (simplify (powerx factor power))))
  287.  
  288.          ;; This does (A+B)^^2 --> A^^2 + A.B + B.A + B^^2
  289.          ;; and (A.B)^^2 --> A.B.A.B
  290.  
  291.          ((and (or (mplusp factor)
  292.                (and (not $dotexptsimp) (mnctimesp factor)))
  293.            (fixnump power)
  294.            (not (greaterp power $expop))
  295.            (plusp power))
  296.           (ncmul factor (ncpower factor (f1- power))))
  297.  
  298.          ;; This does the same thing as above for (A+B)^^(-2)
  299.          ;; and (A.B)^^(-2).  Here the "-" operator does the trick
  300.          ;; for us.
  301.  
  302.          ((and (or (mplusp factor)
  303.                (and (not $dotexptsimp) (mnctimesp factor)))
  304.            (fixnump power)
  305.            (not (greaterp (minus power) $expon))
  306.            (minusp power))
  307.           (ncmul (simpnct-invert factor) (ncpower factor (f1+ power))))
  308.          ((product-with-inner-scalarp factor)
  309.           (let ((p-p (partition-product factor)))
  310.            (mul2 (power (muln (car p-p) t) power)
  311.              (ncpower (cdr p-p) power))))
  312.          ((and $dotassoc (mncexptp factor))
  313.           (ncpower (cadr factor) (mul2 (caddr factor) power)))
  314.          (t (eqtest (list '(mncexpt) factor power) check)))))
  315.  
  316. (declare-top(muzzled nil))
  317.  
  318. (defun simpnct-invert (exp)
  319.        (cond ((mnctimesp exp)
  320.           (ncmuln (nreverse (mapcar #'simpnct-invert (cdr exp))) t))
  321.          ((and (mncexptp exp) (integerp (caddr exp)))
  322.           (ncpower (cadr exp) (minus (caddr exp))))
  323.          (t (list '(mncexpt simp) exp -1))))
  324.  
  325. (defun identitymx (x)
  326.  (if (and ($listp (cadr x)) (= (length (cdr x)) (length (cdadr x))))
  327.      (simplifya (cons (car x) (cdr ($ident (length (cdr x))))) t)
  328.      $dotident))
  329.  
  330. ;  This function incorporates the hairy search which enables such
  331. ;  simplifications as (. a b a b) --> (^^ (. a b) 2).  It assumes
  332. ;  that FIRST-FACTOR is not a dot product and that REMAINDER is.
  333. ;  For the product (. a b c d e), three basic types of comparisons
  334. ;  are done:
  335. ;   
  336. ;  1)  a <---> b        first-factor <---> inner-product
  337. ;      a <---> (. b c)
  338. ;      a <---> (. b c d)
  339. ;      a <---> (. b c d e)    (this case handled in SIMPNCT)
  340. ;   
  341. ;  2) (. a b) <---> c        outer-product <---> (car rest)
  342. ;     (. a b c) <---> d
  343. ;     (. a b c d) <---> e
  344. ;   
  345. ;  3) (. a b) <---> (. c d)    outer-product <---> (firstn rest)
  346. ;   
  347. ;  Note that INNER-PRODUCT and OUTER-PRODUCT share list structure which
  348. ;  is clobbered as new terms are added.
  349.  
  350. (defun simpnct-merge-product (first-factor remainder)
  351.   (let ((half-product-length (// (f1+ (length remainder)) 2))
  352.     (inner-product (car remainder))
  353.     (outer-product (list '(mnctimes) first-factor (car remainder))))
  354.        (do ((merge-length 2 (f1+ merge-length))
  355.         (rest (cdr remainder) (cdr rest)))
  356.        ((null rest) outer-product)
  357.        (cond ((simpnct-alike first-factor inner-product)
  358.           (return
  359.            (ncmuln
  360.             (cons (simpnct-merge-factors first-factor inner-product)
  361.               rest)
  362.             t)))
  363.          ((simpnct-alike outer-product (car rest))
  364.           (return
  365.            (ncmuln
  366.             (cons (simpnct-merge-factors outer-product (car rest))
  367.               (cdr rest))
  368.             t)))
  369.          ((and (not (> merge-length half-product-length))
  370.                (alike1 outer-product
  371.                    (cons '(mnctimes)
  372.                      (firstn merge-length rest))))
  373.           (return
  374.            (ncmuln (cons (ncpower outer-product 2)
  375.                  (nthcdr merge-length rest))
  376.                t)))
  377.          ((= merge-length 2)
  378.           (setq inner-product
  379.             (cons '(mnctimes) (cddr outer-product)))))
  380.        (rplacd (last inner-product) (ncons (car rest))))))
  381.  
  382. (defun simpnct-antisym-check (l check)
  383.   (let (sign)
  384.        (cond ((and (get 'mnctimes '$antisymmetric) (cddr l))
  385.           (setq l (bbsort1 l))
  386.           (cond ((equal l 0) 0)
  387.             ((prog1 (null sign)
  388.                 (setq l (eqtest (cons '(mnctimes) l) check)))
  389.              l)
  390.             (t (neg l))))
  391.          (t (eqtest (cons '(mnctimes) l) check)))))
  392.